home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / blas / ztrsm.f < prev    next >
Text File  |  1996-07-19  |  14KB  |  415 lines
  1.       SUBROUTINE ZTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
  2.      $                   B, LDB )
  3. *     .. Scalar Arguments ..
  4.       CHARACTER*1        SIDE, UPLO, TRANSA, DIAG
  5.       INTEGER            M, N, LDA, LDB
  6.       COMPLEX*16         ALPHA
  7. *     .. Array Arguments ..
  8.       COMPLEX*16         A( LDA, * ), B( LDB, * )
  9. *     ..
  10. *
  11. *  Purpose
  12. *  =======
  13. *
  14. *  ZTRSM  solves one of the matrix equations
  15. *
  16. *     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
  17. *
  18. *  where alpha is a scalar, X and B are m by n matrices, A is a unit, or
  19. *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
  20. *
  21. *     op( A ) = A   or   op( A ) = A'   or   op( A ) = conjg( A' ).
  22. *
  23. *  The matrix X is overwritten on B.
  24. *
  25. *  Parameters
  26. *  ==========
  27. *
  28. *  SIDE   - CHARACTER*1.
  29. *           On entry, SIDE specifies whether op( A ) appears on the left
  30. *           or right of X as follows:
  31. *
  32. *              SIDE = 'L' or 'l'   op( A )*X = alpha*B.
  33. *
  34. *              SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
  35. *
  36. *           Unchanged on exit.
  37. *
  38. *  UPLO   - CHARACTER*1.
  39. *           On entry, UPLO specifies whether the matrix A is an upper or
  40. *           lower triangular matrix as follows:
  41. *
  42. *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
  43. *
  44. *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
  45. *
  46. *           Unchanged on exit.
  47. *
  48. *  TRANSA - CHARACTER*1.
  49. *           On entry, TRANSA specifies the form of op( A ) to be used in
  50. *           the matrix multiplication as follows:
  51. *
  52. *              TRANSA = 'N' or 'n'   op( A ) = A.
  53. *
  54. *              TRANSA = 'T' or 't'   op( A ) = A'.
  55. *
  56. *              TRANSA = 'C' or 'c'   op( A ) = conjg( A' ).
  57. *
  58. *           Unchanged on exit.
  59. *
  60. *  DIAG   - CHARACTER*1.
  61. *           On entry, DIAG specifies whether or not A is unit triangular
  62. *           as follows:
  63. *
  64. *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
  65. *
  66. *              DIAG = 'N' or 'n'   A is not assumed to be unit
  67. *                                  triangular.
  68. *
  69. *           Unchanged on exit.
  70. *
  71. *  M      - INTEGER.
  72. *           On entry, M specifies the number of rows of B. M must be at
  73. *           least zero.
  74. *           Unchanged on exit.
  75. *
  76. *  N      - INTEGER.
  77. *           On entry, N specifies the number of columns of B.  N must be
  78. *           at least zero.
  79. *           Unchanged on exit.
  80. *
  81. *  ALPHA  - COMPLEX*16      .
  82. *           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
  83. *           zero then  A is not referenced and  B need not be set before
  84. *           entry.
  85. *           Unchanged on exit.
  86. *
  87. *  A      - COMPLEX*16       array of DIMENSION ( LDA, k ), where k is m
  88. *           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
  89. *           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
  90. *           upper triangular part of the array  A must contain the upper
  91. *           triangular matrix  and the strictly lower triangular part of
  92. *           A is not referenced.
  93. *           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
  94. *           lower triangular part of the array  A must contain the lower
  95. *           triangular matrix  and the strictly upper triangular part of
  96. *           A is not referenced.
  97. *           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
  98. *           A  are not referenced either,  but are assumed to be  unity.
  99. *           Unchanged on exit.
  100. *
  101. *  LDA    - INTEGER.
  102. *           On entry, LDA specifies the first dimension of A as declared
  103. *           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
  104. *           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
  105. *           then LDA must be at least max( 1, n ).
  106. *           Unchanged on exit.
  107. *
  108. *  B      - COMPLEX*16       array of DIMENSION ( LDB, n ).
  109. *           Before entry,  the leading  m by n part of the array  B must
  110. *           contain  the  right-hand  side  matrix  B,  and  on exit  is
  111. *           overwritten by the solution matrix  X.
  112. *
  113. *  LDB    - INTEGER.
  114. *           On entry, LDB specifies the first dimension of B as declared
  115. *           in  the  calling  (sub)  program.   LDB  must  be  at  least
  116. *           max( 1, m ).
  117. *           Unchanged on exit.
  118. *
  119. *
  120. *  Level 3 Blas routine.
  121. *
  122. *  -- Written on 8-February-1989.
  123. *     Jack Dongarra, Argonne National Laboratory.
  124. *     Iain Duff, AERE Harwell.
  125. *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  126. *     Sven Hammarling, Numerical Algorithms Group Ltd.
  127. *
  128. *
  129. *     .. External Functions ..
  130.       LOGICAL            LSAME
  131.       EXTERNAL           LSAME
  132. *     .. External Subroutines ..
  133.       EXTERNAL           XERBLA
  134. *     .. Intrinsic Functions ..
  135.       INTRINSIC          DCONJG, MAX
  136. *     .. Local Scalars ..
  137.       LOGICAL            LSIDE, NOCONJ, NOUNIT, UPPER
  138.       INTEGER            I, INFO, J, K, NROWA
  139.       COMPLEX*16         TEMP
  140. *     .. Parameters ..
  141.       COMPLEX*16         ONE
  142.       PARAMETER        ( ONE  = ( 1.0D+0, 0.0D+0 ) )
  143.       COMPLEX*16         ZERO
  144.       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  145. *     ..
  146. *     .. Executable Statements ..
  147. *
  148. *     Test the input parameters.
  149. *
  150.       LSIDE  = LSAME( SIDE  , 'L' )
  151.       IF( LSIDE )THEN
  152.          NROWA = M
  153.       ELSE
  154.          NROWA = N
  155.       END IF
  156.       NOCONJ = LSAME( TRANSA, 'T' )
  157.       NOUNIT = LSAME( DIAG  , 'N' )
  158.       UPPER  = LSAME( UPLO  , 'U' )
  159. *
  160.       INFO   = 0
  161.       IF(      ( .NOT.LSIDE                ).AND.
  162.      $         ( .NOT.LSAME( SIDE  , 'R' ) )      )THEN
  163.          INFO = 1
  164.       ELSE IF( ( .NOT.UPPER                ).AND.
  165.      $         ( .NOT.LSAME( UPLO  , 'L' ) )      )THEN
  166.          INFO = 2
  167.       ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
  168.      $         ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
  169.      $         ( .NOT.LSAME( TRANSA, 'C' ) )      )THEN
  170.          INFO = 3
  171.       ELSE IF( ( .NOT.LSAME( DIAG  , 'U' ) ).AND.
  172.      $         ( .NOT.LSAME( DIAG  , 'N' ) )      )THEN
  173.          INFO = 4
  174.       ELSE IF( M  .LT.0               )THEN
  175.          INFO = 5
  176.       ELSE IF( N  .LT.0               )THEN
  177.          INFO = 6
  178.       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
  179.          INFO = 9
  180.       ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
  181.          INFO = 11
  182.       END IF
  183.       IF( INFO.NE.0 )THEN
  184.          CALL XERBLA( 'ZTRSM ', INFO )
  185.          RETURN
  186.       END IF
  187. *
  188. *     Quick return if possible.
  189. *
  190.       IF( N.EQ.0 )
  191.      $   RETURN
  192. *
  193. *     And when  alpha.eq.zero.
  194. *
  195.       IF( ALPHA.EQ.ZERO )THEN
  196.          DO 20, J = 1, N
  197.             DO 10, I = 1, M
  198.                B( I, J ) = ZERO
  199.    10       CONTINUE
  200.    20    CONTINUE
  201.          RETURN
  202.       END IF
  203. *
  204. *     Start the operations.
  205. *
  206.       IF( LSIDE )THEN
  207.          IF( LSAME( TRANSA, 'N' ) )THEN
  208. *
  209. *           Form  B := alpha*inv( A )*B.
  210. *
  211.             IF( UPPER )THEN
  212.                DO 60, J = 1, N
  213.                   IF( ALPHA.NE.ONE )THEN
  214.                      DO 30, I = 1, M
  215.                         B( I, J ) = ALPHA*B( I, J )
  216.    30                CONTINUE
  217.                   END IF
  218.                   DO 50, K = M, 1, -1
  219.                      IF( B( K, J ).NE.ZERO )THEN
  220.                         IF( NOUNIT )
  221.      $                     B( K, J ) = B( K, J )/A( K, K )
  222.                         DO 40, I = 1, K - 1
  223.                            B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
  224.    40                   CONTINUE
  225.                      END IF
  226.    50             CONTINUE
  227.    60          CONTINUE
  228.             ELSE
  229.                DO 100, J = 1, N
  230.                   IF( ALPHA.NE.ONE )THEN
  231.                      DO 70, I = 1, M
  232.                         B( I, J ) = ALPHA*B( I, J )
  233.    70                CONTINUE
  234.                   END IF
  235.                   DO 90 K = 1, M
  236.                      IF( B( K, J ).NE.ZERO )THEN
  237.                         IF( NOUNIT )
  238.      $                     B( K, J ) = B( K, J )/A( K, K )
  239.                         DO 80, I = K + 1, M
  240.                            B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
  241.    80                   CONTINUE
  242.                      END IF
  243.    90             CONTINUE
  244.   100          CONTINUE
  245.             END IF
  246.          ELSE
  247. *
  248. *           Form  B := alpha*inv( A' )*B
  249. *           or    B := alpha*inv( conjg( A' ) )*B.
  250. *
  251.             IF( UPPER )THEN
  252.                DO 140, J = 1, N
  253.                   DO 130, I = 1, M
  254.                      TEMP = ALPHA*B( I, J )
  255.                      IF( NOCONJ )THEN
  256.                         DO 110, K = 1, I - 1
  257.                            TEMP = TEMP - A( K, I )*B( K, J )
  258.   110                   CONTINUE
  259.                         IF( NOUNIT )
  260.      $                     TEMP = TEMP/A( I, I )
  261.                      ELSE
  262.                         DO 120, K = 1, I - 1
  263.                            TEMP = TEMP - DCONJG( A( K, I ) )*B( K, J )
  264.   120                   CONTINUE
  265.                         IF( NOUNIT )
  266.      $                     TEMP = TEMP/DCONJG( A( I, I ) )
  267.                      END IF
  268.                      B( I, J ) = TEMP
  269.   130             CONTINUE
  270.   140          CONTINUE
  271.             ELSE
  272.                DO 180, J = 1, N
  273.                   DO 170, I = M, 1, -1
  274.                      TEMP = ALPHA*B( I, J )
  275.                      IF( NOCONJ )THEN
  276.                         DO 150, K = I + 1, M
  277.                            TEMP = TEMP - A( K, I )*B( K, J )
  278.   150                   CONTINUE
  279.                         IF( NOUNIT )
  280.      $                     TEMP = TEMP/A( I, I )
  281.                      ELSE
  282.                         DO 160, K = I + 1, M
  283.                            TEMP = TEMP - DCONJG( A( K, I ) )*B( K, J )
  284.   160                   CONTINUE
  285.                         IF( NOUNIT )
  286.      $                     TEMP = TEMP/DCONJG( A( I, I ) )
  287.                      END IF
  288.                      B( I, J ) = TEMP
  289.   170             CONTINUE
  290.   180          CONTINUE
  291.             END IF
  292.          END IF
  293.       ELSE
  294.          IF( LSAME( TRANSA, 'N' ) )THEN
  295. *
  296. *           Form  B := alpha*B*inv( A ).
  297. *
  298.             IF( UPPER )THEN
  299.                DO 230, J = 1, N
  300.                   IF( ALPHA.NE.ONE )THEN
  301.                      DO 190, I = 1, M
  302.                         B( I, J ) = ALPHA*B( I, J )
  303.   190                CONTINUE
  304.                   END IF
  305.                   DO 210, K = 1, J - 1
  306.                      IF( A( K, J ).NE.ZERO )THEN
  307.                         DO 200, I = 1, M
  308.                            B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
  309.   200                   CONTINUE
  310.                      END IF
  311.   210             CONTINUE
  312.                   IF( NOUNIT )THEN
  313.                      TEMP = ONE/A( J, J )
  314.                      DO 220, I = 1, M
  315.                         B( I, J ) = TEMP*B( I, J )
  316.   220                CONTINUE
  317.                   END IF
  318.   230          CONTINUE
  319.             ELSE
  320.                DO 280, J = N, 1, -1
  321.                   IF( ALPHA.NE.ONE )THEN
  322.                      DO 240, I = 1, M
  323.                         B( I, J ) = ALPHA*B( I, J )
  324.   240                CONTINUE
  325.                   END IF
  326.                   DO 260, K = J + 1, N
  327.                      IF( A( K, J ).NE.ZERO )THEN
  328.                         DO 250, I = 1, M
  329.                            B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
  330.   250                   CONTINUE
  331.                      END IF
  332.   260             CONTINUE
  333.                   IF( NOUNIT )THEN
  334.                      TEMP = ONE/A( J, J )
  335.                      DO 270, I = 1, M
  336.                        B( I, J ) = TEMP*B( I, J )
  337.   270                CONTINUE
  338.                   END IF
  339.   280          CONTINUE
  340.             END IF
  341.          ELSE
  342. *
  343. *           Form  B := alpha*B*inv( A' )
  344. *           or    B := alpha*B*inv( conjg( A' ) ).
  345. *
  346.             IF( UPPER )THEN
  347.                DO 330, K = N, 1, -1
  348.                   IF( NOUNIT )THEN
  349.                      IF( NOCONJ )THEN
  350.                         TEMP = ONE/A( K, K )
  351.                      ELSE
  352.                         TEMP = ONE/DCONJG( A( K, K ) )
  353.                      END IF
  354.                      DO 290, I = 1, M
  355.                         B( I, K ) = TEMP*B( I, K )
  356.   290                CONTINUE
  357.                   END IF
  358.                   DO 310, J = 1, K - 1
  359.                      IF( A( J, K ).NE.ZERO )THEN
  360.                         IF( NOCONJ )THEN
  361.                            TEMP = A( J, K )
  362.                         ELSE
  363.                            TEMP = DCONJG( A( J, K ) )
  364.                         END IF
  365.                         DO 300, I = 1, M
  366.                            B( I, J ) = B( I, J ) - TEMP*B( I, K )
  367.   300                   CONTINUE
  368.                      END IF
  369.   310             CONTINUE
  370.                   IF( ALPHA.NE.ONE )THEN
  371.                      DO 320, I = 1, M
  372.                         B( I, K ) = ALPHA*B( I, K )
  373.   320                CONTINUE
  374.                   END IF
  375.   330          CONTINUE
  376.             ELSE
  377.                DO 380, K = 1, N
  378.                   IF( NOUNIT )THEN
  379.                      IF( NOCONJ )THEN
  380.                         TEMP = ONE/A( K, K )
  381.                      ELSE
  382.                         TEMP = ONE/DCONJG( A( K, K ) )
  383.                      END IF
  384.                      DO 340, I = 1, M
  385.                         B( I, K ) = TEMP*B( I, K )
  386.   340                CONTINUE
  387.                   END IF
  388.                   DO 360, J = K + 1, N
  389.                      IF( A( J, K ).NE.ZERO )THEN
  390.                         IF( NOCONJ )THEN
  391.                            TEMP = A( J, K )
  392.                         ELSE
  393.                            TEMP = DCONJG( A( J, K ) )
  394.                         END IF
  395.                         DO 350, I = 1, M
  396.                            B( I, J ) = B( I, J ) - TEMP*B( I, K )
  397.   350                   CONTINUE
  398.                      END IF
  399.   360             CONTINUE
  400.                   IF( ALPHA.NE.ONE )THEN
  401.                      DO 370, I = 1, M
  402.                         B( I, K ) = ALPHA*B( I, K )
  403.   370                CONTINUE
  404.                   END IF
  405.   380          CONTINUE
  406.             END IF
  407.          END IF
  408.       END IF
  409. *
  410.       RETURN
  411. *
  412. *     End of ZTRSM .
  413. *
  414.       END
  415.